/*  combD.cpp by K.Tsuru  */ 
/**************************************
It provides a binomial coefficient nCk
 double version. When n < nMax_combD(=48) does not overflow.
See combL(n,k) in "slfcombl.cpp" for argorithm.
**************************************/
#include "tools.h"
#define USE_modf_IN_combD 0

double combD(ulong n, ulong k){
    if(n > nMax_combD) return 0.0;
    if(n < k) return 0.0;
    k = (n < n - k) ? n : n-k; // if n==k, k becomes 0.
    if(k == 0) return 1.0; // includes the case n == k
#if USE_modf_IN_combD
    double u = n;               // work variable in double
    double umax = INT_MAX_DBL;  // = 10^DOUBLE_FIG
    ulong j;
    for(j = 1; j< k; j++){
        u *= n -j;
        if(u >= umax) return 0.0; // overflow
        u /= j+1;
    }
    double ip; // integer part
    if(modf(u, &ip) > 0.0) return 0.0; // not integer
    return ip;
#else
    double u = n;   // work variable in double
    ulong j;
    for(j = 1; j< k; j++){
        u *= n -j;
        u /= j+1;
    }
    return u;
#endif
}